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Abstract. We suggest a two-dimensional wavelet devised to deduce the large-scale structure of a physical field 
(e.g., the Galactic magnetic field) from its integrals along straight paths from irregularly spaced data points to a 
fixed interior point (the observer). The method can be applied to the analysis of pulsar rotation and dispersion 
measures in terms of the large-scale Galactic magnetic field and electron density. The method does not use any d 
priori assumptions about the physical field and can be considered as an algorithm of wavelet differentiation. We 
argue that a certain combination of the wavelet transformation with model fitting would be most efficient in the 
interpretation of the available pulsar RM data. 
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1. Introduction 


Tomography, understood broadly, is a reconstruction of 
a multidimensional physical field from its integral pro¬ 
jections obtained by exposing it in different aspects. 
Tomography problems often occur in astronomy, espe¬ 
cially radio astronomy, where the intervening medium is 
transparent and the observable quantities represent cer¬ 
tain integrals along the line of sight. An important ex¬ 
ample is the studies of the magneto-ionic medium in the 
Milky Way using the Faraday rotation measure RM and 
dispersion measure DM of pulsars - both are integrals 
along the line of sight involving interstellar magnetic field 
and thermal electron density. A peculiar feature of the as¬ 
tronomical tomography is that just one vantage point is 
available, as the observer is confined to the close vicinity 
of the Sun. Astronomical tomography then focuses on ex¬ 
tended objects which can be meaningfully analyzed using 
the integral projections obtained in a variety of directions. 


In this paper we suggest a method to assess the spa¬ 
tial structure of the global Galactic magnetic field using 
pulsar RM (and also DM as the Faraday rotation mea¬ 
sure depends on the thermal electron density), one of the 
most informative tracers of the large-scale component of 
the magnetic field (see e.g., Beck 2000). A fundamental 
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problem here is that the derivation of magnetic field B 
from the observable integral 


RM = A / UeB ■ ds , 

Jo 


( 1 ) 


in fact involves differentiation of the observed quantity 
with respect to r resulting in a catastrophic amplification 
of noise. Moreover, the distances to pulsars, r, are contam¬ 
inates by large errors, both random and systematic. The 
problem is further aggravated by the fact that the pul¬ 
sars (and extragalactic radio sources), that provide the 
RM data, are irregularly spaced and cover the sky very 
inhomogeneously. Here rig is the number density of ther¬ 
mal electrons, r is the distance to the radio source and 
K = 0.81 radm“^ cm^ ^G“^ pc“^; the positive direction 
of s is that towards the observer, so that a field directed 
towards the observer produces positive RM. The thermal 
electron density can be obtained from the dispersion mea¬ 
sure of pulsars ( Manchester 1972 , 1974 ): 


DM = 


Tie ds , 


( 2 ) 


and this involves similar complications. 

The idea of this paper is twofold. Firstly, we use 
wavelets to filter out the noise in the RM data result¬ 
ing from small-scale fluctuations in the magneto-ionic 
medium and from the irregular spacing of the data 
points. Secondly, we introduce a new, specialized family of 
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Fig. 1. The 323 pulsars with known RM from the cat- Fig. 2. The thermal electron density rie distribution in the 
alogue of Faylor et al. (199^ are shown in the galactic Galac tic plane, following the model of Taylor & Cordes 
plane. Here and below, scale is given in kiloparsecs and (1993|), with the Galactic centre marked with cross and 


the Galactic centre is shown with cross; the Sun is at the 
centre of the frame. 


the Sun at the centre of the frame. Darker shades of grey 
correspond to larger densities. 


wavelets devised to avoid (or at least to minimize) noise 
amplification in dealing with an observable represented by 
an integral along the line of sight, given that the available 
data probe a range of distances into the field localization 
region. The Faraday rotation and dispersion measures of 
pulsars are suitable observables for such an analysis. 

The plan of the paper is as follows. In Sect. ||we discuss 
the available samples of pulsar RM and DM and methods 
used to obtain B and rie from them. Section introduces 
the basic ideas of our method, which can be applied in 
various contexts. The algorithm is described in Sect. 
and the advantages and limitations of the method are dis¬ 
cussed in Sect. |^. The efficiency of the method when ap¬ 
plied to the pulsar RM is discussed in Sect. 

2. Obtaining magnetic field from RM and 
electron density from DM 

Faraday rotation measures of pulsars probe the Galactic 
magnetic field in a variety of directions. A pulsar’s own 
contribution to the observed RM is minor because the pul¬ 
sar magnetosphere is populated by electron-positron pairs 
resulti ng in zero net Faraday rotation. This makes pulsars 


netic field of the Mi lky Way (e.g., Rand & Kulkarni 198S; 
Rand fc Lyne 1994 ), especially in conjunction with their 
dispersion measures. However, there are several compli¬ 
cations in the derivation of magnetic field and electron 
density from pulsar RM and DM. 

Firstly, pulsars are very non-uniformly distributed in 
space. They concentrate in spiral arms and the number 
density of known pulsars decreases rapidly with distance 
from the Sun as shown in Fig. ^ Even though the number 
of known pulsars rapidly grows with time, the non-uniform 


nature of their spatial distribution remains. Most pulsars 
are found near the Galactic midplane. Yet, some of them 
are located far away from the midplane, and thereby allow 
a study of the vertical distribution of the magneto-ionic 
medium. Nevertheless, in this paper we restrict ourselves 
to a two-dimensional analysis by projecting the pulsar 
data onto the galactic plane; our method can be general¬ 
ized to three dimensions, but better data would be needed 
to justify this. 

Secondly, the electron density can be obtained from 
dispersion measures only if the distance to the pulsars is 
known. Distance estimates now exist for a few hundreds 
of pulsars, resulting from three basic techniques: neutral 
hydrogen absorption (in combination with the Galactic 
rotation curve), trigonometric parallax and from associa¬ 
tions with objects of known distance ( Lorimer 2001 ). The 
distance of a pulsar can b e obtained from DM if the dis¬ 
tribution of rie is known. Taylor & Cordes (1993), based 
mainly on the scintillation and dispersion measure data, 
have suggested an electron density model expected to pro¬ 
vide distance estimates with an average uncertainty of 
about 30%. However, distances to individual pulsars may 


be as wrong as by as a factor of two (e.g., lohnston et al 


a majcr source of information about the large-scale mag- 2001). 


The method suggested here can be applied to pulsar 
dispersion measures to obtain thermal electron density us¬ 
ing pulsars with known distances. Then magnetic held can 
be obtained from RM using the same method. However, 
dispersion measures alone are not sufficient to obtain a 
reliable distribution of rig, and a model for the electron 
density should be used in order to derive magnetic held 
distribution from Faraday rotation measures. 

The largest catalogue of pulsars contains 706 objects 


(Taylor et al. 1995). Faraday rotation and dispersion mea- 
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sures, together with the distances are known for only 323 
of them. 

The simplest estimate of the line of sight component of 
the large-scale magnetic field B can be obtained by divid¬ 
ing RM by DM. This would be justified if both B and rie 
were uniform along the line of sight, which is definitely not 
the case. A more elaborate analysis involves a pair of pul¬ 
sars located close to each other in the sky, so Ue between 
the two pulsars can be obtained from the increments in 
DM and r as Ue — ADM/Ar. Then the longitudinal mag¬ 
netic field follows as i3|| ~ ARM/(ArAne). However, Ar 
is most often too large to make the resulting estimate re¬ 
liable. The nature of the problem is intrinsically related 
to that arising in differentiating unevenly sampled data 
(cf. Ruzmaikin et al. 1988, p. 34). A possible resolution 
is to use multivariate statistical analysis where a model 
for magnetic held is adopted and its parameters are ob¬ 


tained by htting (Ruzmaikin et al. 1988; Rand & Kulkarni 
1989|^^nd fc Lyne 1994|). A disadvantage of this (and any 


other) analysis of this kind is that one has to impose d 
priori model restrictions on the structure of the magnetic 
held. Furthermore, it is often difficult, if not impossible, 
to get convincing evidence that the model adopted is com¬ 
patible with observations because errors in the observable 
quantities are rarely known conhdently and this hampers 
the application of statistical tests for the quality of the 
ht. Here we attempt at developing a method of held re¬ 
construction from the integrals ([I]) and (H) which allows 
us to avoid these shortcomings. Our method is based on 
wavelet analysis. 


3. The wavelet analysis 

The data similar to those presented in Fig. [| obviously 
would not allow one to reconstruct the held in all de¬ 
tails. What is possible, however, is to reveal structures 
whose scale is larger than the size of the gaps between 
the points observed. For this purpose, one should sepa¬ 
rate large-scale structures from noise resulting not only 
from statistical errors but also from the irregular data 
grid. Recently wavelets have become an efficient tool for 
signal analysis, especially useful in scale separation prob¬ 
lems. The main idea of this approach is a decomposi¬ 
tion of the signal over a basis of self-similar functions, 
known as wavelets, that are localized in both physical and 
Fourier space. The wavelet transform can be considered 
as a generalization of the Fourier transform, which allows 
the derivation of the local spectral properties of the sig¬ 
nal. Comprehensive reviews of the wavelet analysis can be 
found in Farge et al. (1993 ); Grossmann fc Morlet (1984 ); 
Holschneider (1995 ). 

In terms of the polar coordinates (r,ip), the continu¬ 
ous two-dimensional wavelet transform (also known as the 
wavelet coefficient) of a function F(r, tp) is given by 


TT OO 

'j{a,r,ip) = ^ j j 'tjja,r,v{r',ip')F{r',ip')r'dr'd(p', 

— 7T 0 


( 3 ) 


where 'ipa,r,ipi'r', is a family of wavelets characterized 
by position in the physical space (r,ip) and the scale a 
which can be identified with the radius of the region where 
,(p') = 0(1). A self-similar family of wavelets is 
generated from an ‘analyzing wavelet’ 'I'(r, tp) by transla¬ 
tion and dilation. Many various wavelets have been sug¬ 
gested. The success of wavelet analysis depends on ap¬ 
propriate choice of the wavelet 'I'(r, ip). Most importantly, 
some wavelets provide better resolution if the wave num¬ 
ber (Fourier) space (at the expense of the localization 
information), whereas others are advantageous if higher 
spatial resolution is needed (but then the scale separa¬ 
tion power can be only modest). The wavelet used here 
derives from an isotropic wavelet known as the ‘Mexican 
hat’, which has good resolution in the physical space, 

4'(x) = (2-a:^)exp(-ia;2) , (4) 


where x is the distance from the wavelet center. An 
isotropic wavelet has no preferred direction, and so ro¬ 
tation is not included into the transformations generating 
the wavelet family, which is defined as 




s(r, ip, r',ip') 


( 5 ) 


where s(r, ip, r', ip') is the distance between the (r, ip) (the 
centre of the wavelet) and the current point (r', ip'). 

Ideally, wavelet analysis consists of decomposing the 
signal into contributions from a range of scales, filtering 
out the smaller scales dominated by noise and the recon¬ 
struction of the filtered signal (presumably free of noise) 
using the inverse wavelet transform. However, this pro¬ 
cedure is rarely followed completely because this would 
require observational data of unrealistically high qual¬ 
ity. Instead, we shall consider the wavelet coefficients 
themselves, but at a scale distinguished by the data. 
Specihcally, we focus on those scale(s) which make a dom¬ 
inant contribution to the ‘spectral energy density’ (see 
below). The wavelet coefficients form a three-parametric 
family. For representation and interpretation purposes, we 
shall usually consider a set of two-dimensional maps of 
w[a, r, ip) for a hxed a. 

A useful integral diagnostic of the importance of dif¬ 
ferent scales in the data is the wavelet spectrum, defined 
as 

TT OO 

M{a) = j j w‘^{a,r,ip)r dr dip . (6) 

— TT 0 


In the Fourier analysis, a similar quantity is known as 
the spectral energy density. The normalization factor in 
Eq. (|^ is chosen so that structures of the same amplitude 
(but different scales) contribute equally to M{a). 


4. Wavelet difrerentiation 

In this section we suggest a general method to analyze 
an observable represented by an integral, with variable 
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limit(s), of the physical field whose reconstruction is the 
goal of the analysis. Deriving magnetic field from pulsar 
RM is a perfect example of such a problem, however the 
method can be used in many other problems of image 
analysis. The main idea of the method is to avoid the 
differentiation of the discrete, noisy data. Instead, we dif¬ 
ferentiate the wavelet itself, which is given in an analytical 
form. Thus, the algorithm can be considered as a sort of 
wavelet interpolation used to derive the derivative of an 
observable quantity. 


4.1. The differentiating wavelet 

Let F{r, (p) be an unknown function whose integral 


r 

G{r,ip) = j F{r',ip)dr , 


(7) 


is known for various values of r. In the context of this 
paper, G{r,(p) can be DM or RM (and F{r,(p) is related 
to rig or ngi^ii, respectively), and the origin of the reference 
frame is at the solar position. 

The exact expression for F involves a derivative of G, 
F{r, (fi) = dG{r, (p)jdr. Since G is known only at discrete 
data points, one has to use a smooth interpolation of G 
in order to derive F{r,ip). Random and systematic noise 
makes this procedure unstable. This difficulty also arises in 
standard tomography problems where a two-dimensional 
function is reconstru cted from its integral proje ctions con¬ 
taminated by noise ( Patrickeyev fc Frick 199^ ). 

Wavelet analysis can help to alleviate this problem. 
Given Eq. (0), one can rewrite Eq. (||) as 



w{a, r,ip) = Jil}a,r,v{r', r' dr' dip'. ( 8 ) 

— TT 0 

Integrate by parts to obtain 

TT CX> 

w{a,r,ip) = -^ JJ ,ip')G{r',ip')r'dr'dip' , (9) 

— TT 0 


Fig. 3. (a) The ‘Mexican hat,’ Eq. M, and (b) the as¬ 
sociated differentiating wavelet, Eq. (M). In this illustra¬ 
tion, the wavelets are centered at {r,ip) = (5,0) and have 
the scale a = 2. The zero-level contours are shown solid 
and the dashed contours corresponds to the magnitude of 
1% of the maximum. Darker shades of grey correspond to 
larger values. 


where t') is a new wavelet family related to the 

original wavelet by 


~r ( ' dtpa,r,(p{'r t'P) , <P ) 

, (^ ) = - — - + - - - 


( 10 ) 


Here we have taken into account that G(0, (/?) = 0 and 
that the wavelet is localized in the physical space, 


lim r'ipa,r,v{r',ip') = 0 . (11) 

r'—>oc 

Altogether, we have introduced a new family of 
wavelets that allows us to reconstruct F(r,ip) without 
differentiating the data. The forms of the ‘Mexican hat’ 
and the associated differentiating wavelet are illustrated 
in Fig. H Unlike the ‘Mexican hat’, the differentiating 
wavelet is not isotropic, and has two zeros. 


The new wavelet if has zero mean value, 

TT OO 

j j i’a,T,v(.r',ip')r'dr'dip'= 0 , (12) 

— TT 0 

and so satisfies what is known as the admissibility condi¬ 
tion. Another admissibility condition is that 

TT OO 

J J ip') r'2 dr'dip' = 0 . (13) 

— TT 0 

Equations and imply that the wavelet coefficients 
are independent of a constant component of G(r, ip) and 
its linear trend, respectively, - a behavior required for any 
wavelet in order to have an inverse transformation. Thus, 
F (r, ip) can be determined up to a constant background 
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value proportional to the value of at r = r„iax, where 
’’max is the distance to the most remote pulsar. If r^ax is 
large enough, this background value is close to zero. 

4.2. Implementation for irregularly spaced data 

The calculation of the integral on a discrete, irregu¬ 
lar data grid is not straightforward. The need for a robust 
numerical algorithm is especially demanding since the sep¬ 
aration of pulsars in the available catalogues is often com¬ 
parable to or exceeds the physically interesting scales of 
the large-scale Galactic magnetic field. The difficulties be¬ 
comes crucial if the separation of the data points can be 
comparable to the scale of the wavelet, a. 

In our case, a convenient method to calculate the in¬ 
tegrals is the Monte-Carlo technique, leading to 



-10 -5 0 5 10 


w{a,r,ip) = - 


A5(a, r, (p) 


N 

E 




where N is the number of data points for G(r, ip) and AS 
is the average area element of an individual data point. 
For a uniform distribution of the data points, is a 
constant, but this is not true for a non-homogeneous dis¬ 
tribution. We chose (a function of position and the 
scale of the wavelet, identified with the effective weight of 
data points) as to minimize the error arising from the non- 
uniform distribution of the data points. The contribution 
of a data point to the integral sum ( 0 ) must be small if 
the point is far away from the centre of the wavelet, r' , as 
compared with the scale of the wavelet, a, so that we take 


AS {a, r, (fi) = 


■ N 

E^ 


s{r,(p,ri,ip^] 


(15) 


The value of / AS gives the approximate number of 
sample points within the wavelet localization area. So, 
a'^/AS should be large enough in order to achieve suf¬ 
ficient accuracy. 

Another problem arising in the numerical realization 
of the wavelet transformation is the necessity to satisfy 
the admissibility conditions (0) and (0) with sufficient 
accuracy. We use for this purpose the gapp ed (or adap¬ 
tive) wavelet technique (see Frick et al. 1998| ) wherein the 
‘analyzing wavelet’ is represented in the form 

'I'(x) = h{x)^{x) , (16) 

where <l>(x) is a positive definite function that deter¬ 
mines the position of the wavelet (the envelope). For the 
‘Mexican hat’, Eq. (||), $(a;) = exp(—^a;^) and h{x) = 
2 — x'^. The main idea of the technique is to modify h{x) 
in order to satisfy the two conditions (12) and (p^. This 
can be achieved by introducing two quantities Ci and C 2 , 
functions of scale a and position (r, ip), such that 

'4){x) = [h{x) + Co + Cix\^{x) . (17) 

These functions are determined from the discrete versions 
of Eqs. (12) and for each scale a and wavelet position 
(a V?)- 
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Fig. 4. (a) The test function F and (b) the corresponding 
distribution of G calculated from F using Eq. With 
F = Tie, G is proportional to DM. Lighter shades of grey 
correspond to smaller values. 


5. Testing the method 

In order to assess the possibilities of the method, we con¬ 
sider a test function F{r,ip) = rie similar to that of the 
Milky Way Taylor fc Cordes (1993| ) . We first calculate the 
corresponding values of G{r, tp) cx DM(r, ip) using Eq. (|^ 
and then apply our method to obtain F{r, ip) from G(r, ip) 
on a selection of data grids. 

We show in the upper panel of Fig. ^ a test function 
F = Tie that includes an annular and spiral structures with 
characteristic scales (transverse half-widths) a = 3kpc 
and 0 = 1 kpc, respectively. This distribution in similar 
to (but not identical with) that of Fig. The lower panel 
of Fig. H shows the corresponding G = DM(r, ip) calcu¬ 
lated in the whole plane. 

The quality of the reconstruction of F from G obvi¬ 
ously depends on how well G is sampled. For a reliable 
reconstruction of a simple, isotropic isolated object whose 
structure is similar to that of the wavelet itself, one needs 
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Fig. 5. The wavelet transform of the distribution of 
Fig. ^ on a regular 40 x 40 grid at a scale (a) a = 3kpc 
(the scale of the annulus) and (b) a = 1 kpc (the scale of 
the spirals). 


at least 10 data points within it. Then one would need 
about 2000 data points distributed uniformly in (r, cp) in 
order to detect structures of the smallest scale of Ikpc 
in the 20 kpc x 20 kpc Solar vicinity. However, Faraday 
rotation measures are only known for about 300 pulsars; 
furthermore, their distribution is very non-uniform being 
very sparse beyond 3-5 kpc from the Sun, especially in the 
fourth Galactic quadrant. Thus, the available pulsar data 
can only yield reliable information about structures in a 
narrower vicinity of the Sun and/or at a scale of a few 
kiloparsecs. 

Now we consider a series of wavelet transforms on dif¬ 
ferent data grids, from an ‘ideal’ to a realistic case through 
an ‘acceptable’ one. 

A regular grid is the best for the calculation of the in¬ 
tegral (^). The data point separation of 0.5kpc (needed 
to recover structures at a scale of 1 kpc) requires a sample 
of 1600 points. The wavelet transforms of the distribu¬ 
tion of Fig. ^ sampled on a regular grid of 1600 points 



Fig. 6. The integral wavelet spectrum M{a) for the dis¬ 
tribution of Fig. |b, sampled on a regular 40 x 40 grid, 
with a in kpc. 


are shown in Fig. ^ for the scales a = 3kpc and Ikpc. 
As expected, the annular structure is better pronounced 
at a = 3 kpc (Fig. |a) and the spirals are better visible at 
a = Ikpc (Fig. ||b). Figures ||a and b together reproduce 
the distribution of Fig. ^ reasonably well, demonstrating 
the possibilities of the wavelet transform in scale separa¬ 
tion. However, the scale resolution of the wavelet, derived 
from the ‘Mexican hat’, is only modest, and this contam¬ 
inates the wavelet transform at a = 3 kpc with a trace of 
the spiral arms, and that at a = 1 kpc with a weak signa¬ 
ture of the annulus. The magnitude of the wavelet coeffi¬ 
cient is maximum when the scale of the wavelet is equal to 
the scale of the structure. The integral wavelet spectrum 
M (a), defined in Eq. (H) and shown in Fig. H confirms that 
the scales 1 < a < 3 are dominant in G(r, i^), but they are 
not separated from each other. 


Further, we consider the data of Fig. ^ sampled at the 
same number of points (1600), but now the data points 
are scattered randomly and statistically uniformly in the 
same region 20 kpc x 20 kpc around the Sun. The distor¬ 
tions in the resulting wavelet transform maps of Fig. 
are significant even though the two basic components of 
the distribution are clearly recognizable. The structures in 
the wavelet transform become patchy, especially at small 
scales, because of the large gaps in the data grid. 

Finally, the wa velet transform on the data grid of the 
pulsar catalogue of Taylor et al. (1995| ) (using only the po¬ 
sitions of 323 pulsars with known RM) is shown in Fig. |[ 
The data points are crowded within 3 kpc from the Sun. 
We consider the wavelet transform to be unreliable where 
l^S > a^/5; these regions are marked with uniform grey 
shade in Fig. The spiral arms segments of the original 
distribution of Fig. ^ is recognizable in the 3 kpc vicinity 
of the Sun, but they appear patchy and discontinuous. The 
annulus at the galactocentric distance of 4 kpc is hardly 
visible at all, and the arm-interarm contrast is overesti¬ 
mated. 
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Fig. 7. The wavelet transforms of the distribution of 
Fig. sampled at 1600 randomly distributed points at 
the scales (a) a = 3 kpc and (b) a = 1 kpc. 


6. Discussion 


We have introduced a new wavelet devised for the anal¬ 
ysis of the pulsar Faraday rotation measures in terms of 
the large-scale magnetic field (or of any other observable 
that is an integral of the quantity studied, e.g., the pul¬ 
sar dispersion measures). This is a tomography approach 
because the field is directly reconstructed from its inte¬ 
gral estimator. The method works well with data given on 
a regular mesh or scattered randomly but with gaps be¬ 
tween the data points not exceeding a/2, with a the scale 
of the wavelet. The separation of pulsars with known RM 
exceeds this limit beyond about 3 kpc from the Sun. 

An advantage of the method is that it involves nei¬ 
ther ad hoc assumptions about magnetic field structure 
nor model fitting. However, the wavelet transforms on the 
data grid of the pulsar catalogues available appear rather 
confusing and difficult to interpret. The advantages of 
the wavelet analysis and model fitting can be combined 


in a single approach applied by Frick et al. (2001) to 



Fig. 8. The wavelet transform of the distribution of 
Fig. |b on the real, irregular data grid of the pulsar 
catalogue of Taylor et al. (1995 ) (323 points) at a scale 
a = 1.5 kpc. The region where AS” > a^/5 and the results 
are inaccurate is filled with uniform grey. 


the Faraday rotation measures of extragalactic sources. 
Instead of htting a magnetic field model to noisy RM data, 
these authors fitted the wavelet transform of the model 
RM to the wavelet transform of the observed RM with 
smaller scales (mainly responsible for the nose) filtered 
out. This has resulted in a significant improvement in the 
quality of the analysis. The application of the wavelet in¬ 
troduced here will allow us to fit the wavelet transforms 
of the magnetic field derived from RM to the model. We 
expect that will result in a significant improvement of the 
results against the usual procedure of fitting model RM 
to the observed noisy data. 

Another way to improve the results would be to ana¬ 
lyze simultaneously the Faraday rotation measures of pul¬ 
sars and extragalactic radio sources. Since many extra- 
galactic sources occur at high Galactic latitudes, a plau¬ 
sible assumption about the vertical distribution of the 
magneto-ionic medium will have to be adopted. 
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